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Abstract 

The plane wave stability properties of the conservative schemes of 
Besse and Fei et al. for the cubic Schrodinger equation are analysed. Al- 
though the two methods possess many of the same conservation proper- 
ties, we show that their stability behaviour is very different. An energy 
preserving generalisation of the Fei method with improved stability is 
presented. 



1 Introduction 
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Figure 1.1: The plot shows the numerical solution of the cubic Schrodinger 
equation produced by the method of Fei et al. [5] and the method of Besse [3] , 
in two consecutive time steps. Here uo(x) — e smx i t — 1.9, time step r = 0.01, 
space step h = 27r/1024, A = 2. 



Figure |1.1| shows the numerical approximation in two consecutive time steps 
with two well-known schemes applied to the cubic Schrodinger equation (CSE) . 

\u t + Au = X\u\ 2 u, Aet. (1.1) 

Both schemes have order 2, are symmetric, and conserve mass and energy in a 
way to be made precise below. They are both linearly implicit. Whereas the 



1 



method of Besse exhibits a continuous behaviour, the method of Fei et al. has an 
unstable spurious solution with eigenvalue near —1 which causes the numerical 
solution to alternate in each time step. 

We begin by briefly describing some properties of the CSE and some nu- 
merical methods which have been proposed for approximating its solution. In 



one space dimension, d = 1, (1.1) is completely integrable and for the sake of 
notational simplicity we shall assume in the rest of this paper that there is just 
one space dimension. All the numerical methods we consider can be easily gen- 
eralised to arbitrary d > 1. We may also add that the integrability which is 
particular to the case d = 1 will not be important for the issues discussed here. 
Generally, the density 

p[u] = f \u\ 2 dx, (1.2) 



is a conserved quantity, and it is also usual to work with the Hamiltonian struc- 
ture 

H[u] = J Q|V U | 2 + ^M 4 ) Ax. (1.3) 
In the literature, there exists a large variety of numerical schemes for the inte- 



gration of ( 1.1 ), see for instance [H [H O O (H [H [H [TUl EU - Frequently one sees 



examples of numerical schemes that are conservative, by this we mean that they 



exactly preserve some discretised version of (1.2), ( |l.3[ > , or both. In addition, 
it is often desirable to have schemes that are symmetric, see e.g. Hairer et al. 
[7j for a proper definition. One favourable consequence of having a conservative 
scheme, is that it can be used to control the growth of the numerical solution 
over long times. Generally, most conservative schemes are implicit in the time 
step. But the class of implicit schemes can be further divided into schemes that 
arc fully implicit, and those that are linearly implicit or semi-explicit. In the 
former system of nonlinear equations must be solved in every time step, 

the size of which is equal to the number of spatial degrees of freedom. It is often 
necessary to use a Newton-type iteration for solving this nonlinear system, and 
especially for large time steps the convergence may be slow or the iteration may 
not converge at all. However, using a linearly implicit scheme guarantees that 
the cost is roughly the same in each time step, therefore such schemes are at- 
tractive from the point of view of computational efficiency. Examples of linearly 
implicit schemes which are symmetric and conserve some discretised version of 
the energy when applied to the cubic Schrodinger equation are the method of 
Besse [3] and that of Fei et al. [6] , the latter being derived also in [9] . The former 
scheme is formulated continuously in space, i.e. for each time = to < ^1 < • " " 
the method produces an approximation U n to u{x,t n ) as follows 

An+1/2 , An-1/2 

— = \U n \ 2 , (1.4) 
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u^_-un + A (u^ + u«\ = xr+1/2 fu^ + u . (L5) 



This method can alternatively be written as a two-step scheme, and thus an 
auxiliary approximation is needed for the first step. Besse proposes to set 
([/°, 0" 1 / 2 ) = (u(x, 0), \u(x, 0)| 2 ), and shows that the following two versions 
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of (1.21 and (1.3 1 



D n 



U n \ 2 Ax 7 



H n = J Q|V[/| 2 + ^ n+1/ V l ~ 1/2 ) di, 

are preserved in the sense that D n = D° and E n — E° for all n. 

The method of Fei et al. is different from the Besse method in that it is 
discretised both in time and space, in [B] the scheme is given for one space 
dimension as follows: 



.U, 



n+l 



U" 



2t 



1 /U n+1 4- 77™ _1 
h 2 \ 2 



Alter 



u. 



n+l 



, (1.6) 



jjn 



where 5 X is the centered difference in space, that is, S X U^ = t^™ + i 

Here is an approximation to u(x mi t n ). The associated discretised density 
and energy of this scheme are 

h 



Tjn+1 
u m+l 



- u n+1 


+ 


Jjn _ Jjn 
u m+l u m 


h 




h 



+ x\u, 



" +1 ' 2 |C/"| 2 



In a similar way as for the Besse scheme, (1.6) can be formulated continuously 
in space as 



1 2t 



jjn+l + jjr. 



(1.7) 



The rest of the paper is organised as follows; in sectio n [2] w e analyse the 
plane wave stability properties of the two schemes ( 1.5 ) and ( 1.6 ). We illustrate 



the results using various plots. In section[3]we introduce a two parameter energy 
conserving generalisation of the Fei method that improves stability. 



2 Plane wave solutions, dispersion relations and 
stability 

Most of this section is inspired by the analysis of [TT] where dispersion relations 
for the exact and some numerical solutions of ( 1.1 ) are derived and their linear 
stability is analysed. 



2.1 The exact solution 



The cubic Schrodinger equation (1.1) with periodic boundary conditions sup- 



ports plane wave solutions of the form 

v(x,t) = Q j( kx -» t \ 
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where u> is determined by the dispersion relation 

lu = k 2 + X\a\ 2 . 



(2.1) 



We then consider small perturbations of such a solution, substituting u(x, t) = 
(1 + e(x,t))v(x,t) into (1.1 1, and ignoring 0(e 2 ) terms. We get 

\e t + 2ike x + e xx - X\a\ 2 {e + e) = 0. (2.2) 

The perturbation e{x, t) is periodic, thus we invoke its Fourier expansion 

e(x,t) = Y,£e(t)e Ux . 



Substitute this expansion into (2.2) to obtain 
d 



dt 



e t = i ({-2k£ -I 2 - \\a\ 2 )ei - \\a\ 2 i- t ) 



We take the complex conjugate of this last equation and replace £ by —£, the 
result is the linear system of ODEs 



d_ 
dt 

The eigenvalues of Gi are 







' it ' 


, Gi — 






F-i. 





-£ 2 -2k£- \\a\ 2 -X\a\ 2 

X\a\ 2 £ 2 ~2k£ + X\a\ 2 



X e = (-2k ± ^j£ 2 + 2A|a| 2 ) I 



If A^ is complex, then one eigenvalue of (i Gi) will have a positive real part, and 
the corresponding mode will be unstable. This happens if 

£ 2 < -2A|a| 2 , 

thus instability may only occur when A < 0, usually referred to as the focusing 
case. 



2.2 The scheme of Fei et al. 

We now consider again the scheme 

Tjn+l _ Tjn-1 i /r/™ +1 4- TI n ^ 1 \ o f TJ n+1 A- TJ n ~ 1 \ 

i m 2r m 2 ) =A l^i ( — 2 )■ (2 - 3) 

Substituting a sequence of the form V£ — ae 1 ^ kXm ~ ultnS> in which x m = rah, 
t n = tit, we get the dispersion relation 

tanwr = Xr\a\ 2 + 4psin 2 ~, P=t^- (2.4) 

2 h z 



We perturb this plane wave solution, substituting {/" = (1 + eJJVJJ into (2.3) 
and after ignoring quadratic and higher order terms in we get 

aie m+1 +a e m +a- 1 e m _ 1 - b{e m + e m )~a^ 1 e m+1 ~a e m —ai£ m _ v (Z.b) 
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Here 



ai = pS kh ~" T) 
a = e^ T (i- 2p- q) 
a_i = pe-^ kh+UT \ 



b = 2q cos lot 
a = Arlal 2 



We expand in a Fourier series e™ n = ^£™e lfem and substitute this into (2.5) 
to obtain 

Q e™ +1 = + 1"^,) - Qe"~\ q = aic i£ ' 1 + a + a_i e" itt . 
We now take the complex conjugate of this equation and replace I by — £ to 

i T 



obtain a system of difference equations for E n = 



Eg , £_ 



~ci " 


E n+1 = 


b 


b 


E n - 


Q " 


c- e _ 


b 


b 


c- e 



We find that this difference equation is stable if and only if the polynomial 

p(z) = CfC-tz 4, - b(ct + c^i)z 3 + {c t c_i + c e c_ e )z 2 - b{c_ t + cg)z + c^ t c t (2.6) 

has all its roots in the closed unit disc. Note that p(z) is self-reciprocal, meaning 
that its set of roots is invariant under the transformation z t— > f /z. Each root on 
the unit circle is invariant under this transformation, but any root in the open 
unit disc is mapped to a root outside the unit circle. Thus, for self-reciprocal 
polynomials, stability is equivalent to all roots lying on the unit circle. Note 



that when we use the scheme which is continuous in space (1.7 1, we get again 



the stability polynomial (|2.6|), but where eg is replaced by 



{.-{K + Lf-q) 



K 



The case k = is particularly simple, since the stability polynomial in this case 
will have real coefficients. In that case one may easily derive that p(z) has all 
its root on the unit circle for all values of I if and only if 



A|a| 2 < 



1 — cos h 



Note that this condition is independent of r, and in the limit when h tends to 
zero, one simply gets the condition A|a| 2 < |. This does however not imply 
that the scheme does not converge on finite time intervals [0, t*]. Suppose that 
A|a| 2 > |. For small perturbations, one may expect that the error is roughly 
amplified with a factor in each step that equals the magnitude of the largest 



root of (2.6|, which for k — and in the limit case h — > is of the form 

z* = - (l + rV2AM 2 -^ 2 ) + 0(t 2 ), 

thus locally the size of the spurious solution grows exponentially, the growth 
factor of the unit frequency over the interval [0, t*] being approximately exp(Ci*) 
with C = v^AH 2 - 1. 
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2.3 The scheme of Besse 



In one space dimension the scheme of Besse (1.4), (1.5) is 



-(u 

T 



n+1 



U 



n+1 



u r 



\u n 

\(j) r ' 



u 



n+1 



u n 



(2.7) 
(2.8) 



We use the plane wave solution V n = ae^ kx wt »), which now is continuous in 
space, to get the following relation for to. 



tan — = - [\T\a\ 



As for the Fei scheme we consider the perturbations 

jjn = yn ^ + g nj ( ^n+\ = | fl |2 ^ + § n+^\ 



(2.9) 



Substituting these expressions into (2.7) and ignoring higher order terms yields 
<5™+3 = -<J"-i +2{e n + e n ). (2.10) 



We now plug the last three expressions into (2.8) to obtain 

i ^-iur^+1 _ £ nj + ^-iur f^n+1 + 2 ike^ +1 - (k 2 + A |fl| 2 ) E n+1 

+ \ (e n xx + 2ife< - (k 2 + A |a| 2 ) e n ) = ~A \a\ 2 (e"^ + l) 5 n +i . (2.11) 



Expand e n (x) and 8 n+ ^ (x) in a Fourier series and insert into (2.10) and (2.11 ) 
to get the system 








-b 


" 


















C-t 





-6 


E n+1 = 





-C-i 














1 





2 


2 


-1 














1 




2 


2 





-1 



E n E n 



—n 



where we have defined 



c t = (2i — (L + Kf - q)e~ luJT 
b = 2gcos — , 

and where L = y/r£ and K = y/rk as before. Notice that cp and b are defined 
differently than in the Fei case. We find that the characteristic polynomial is 
(z + l)p(z) with 

p(z) = c e c^ e z 3 + {cic-t + cic-t + c£c-i - 26 (q + C-t)) z 2 

+ (ciC-t + C£C-£ + cic-i - 26 (ce + c-t)) z + c e c^ e . 
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To shorten the notation we divide by cgC-i and define / and g such that 

p(z) = z 3 + gz 2 + gz + f, 

and p(z) has the same roots as p(z). The polynomial is self- reciprocal and we 
have stability only if all three roots lie on the unit circle. We proceed to express 
the stability region S in terms of g for a given /. The key observation is that 
for points on the boundary dS at least two of the roots are coalescing. For such 
values of / and <?, we can write the polynomial as 

By comparing coefficients we get that for a given value of / we can parametrise 
the stability boundary for g as follows 

g(9)=c 2W -2fc- w . 

Since |/| = 1 we have that \g\ < 3 is necessary for stability while \g\ < 1 is 
sufficient. For k = the latter condition becomes £ 2 < — 2 A| a| 2 , which is the 
same as for the exact solution. 

For k = the two roots of p(z) with the largest magnitude are on the form 

z* = 1 ± r<V-2A|a| 2 -£ 2 + 0(t 2 ). 

In this case the eigenvalue is near 1, not —1 as in the Fei case, which explains 
why the Besse scheme does not exhibit an alternating behaviour for k = 0. 



Extension to more space dimensions. The scheme of Fei (2.3) has an 
obvious generalisation to d space dimensions with corresponding energy and 
density functions. We may also consider the general form of the Besse scheme 



(1.4), (1.5), and introduce plane wave solutions in d dimensions 

it(se,t) = ae 1 ^" - "*) 

for space variables x = (xi, . . . , Xd) and wave numbers fc = (kx, . . . , fed). Exact 
solutions of this form satisfy the dispersion relation 

u: = \k\ 2 + \\a\ 2 , 

and these solutions are stable with respect to perturbations 

ee(x,t) = e e (t)e ie - x , 

whenever \£\ 2 > -2A|a| 2 . 

For the method of Fei et al. one finds again that for k = the scheme 
is stable to perturbations only if A|a| 2 < \ in the limit h — > 0, the critical 
perturbation wavenumber vectors being the canonical unit vectors in M. d . The 
dispersion relation is now 



tan W r = Ar|a| 2 +4p^sin 2 ^ ^ Ar|a| 2 + r|fe| J 

3=1 
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Figure 2.1: A comparison of the dispersion relations. 



Also the method of Besse generalizes similarly in more space dimensions, the 
dispersion relation being 

and the stability polynomial is obtained simply by replacing q by 
ct = (2i- r\£ + k\ 2 -dq)e- iuiT . 



2.4 Numerical results 

We compare the schemes of Besse and Fei et al. in the limit h — > 0, in which 



case p(z) depends on the three parameters K, L and q. In figure 2.1 

2 



we compare 



the exact dispersion relation (2.1 1 with Fei (2.4 1 and Besse (2.9). Both relations 
express lot in terms of k 2 + X\a\ 2 and for the Besse method it is obtained by 
replacing the time step r by r/2 in the Fei method. 

Figure |2.2| shows a numerical computation of the stability region for both 
Besse and Fei when fixing K = k — in p(z). These plots are obtained 
from the explicit expressions in the analysis of the previous sections. The case 



K = ky/r = 1 is shown in figure 2.3 For reference, we have included the stabil- 
ity boundary of the exact solution as a broken line. The stability of the Besse 
scheme differs slightly from that of the exact solution when A is negative, how- 
ever, the stability is retained for positive A. The Fei scheme is again unstable 
in the defocusing case. We have observed that when K is further increased, 
then also the Besse method becomes unstable for small modes in the defocusing 
regime. 

For a numerical plane wave solution to be stable, none of its Fourier modes 
I should be amplified by the method. The plot in figure [2~4| shows the stability 
region in the qK-p\&ne for the Besse method. A pair (q, K) is unstable if the 
largest root of the stability polynomial exceeds 1 in modulus for some real value 
of L. It is then only of interest to consider the defocusing case (q > 0) since the 
exact solution itself is unstable in this sense for all negative q. The Fei scheme 
is unstable for all (q, K) where q ^ 0. 
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Besse (K = 0) Fei (K = 0) 
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Figure 2.2: Stability for K = (grey is unstable). 

Besse (K = 1) Fei (K = 1) 
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Figure 2.3: Stability for K = 1. 



Besse 




Figure 2.4: Stability for all I. 
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3 An energy preserving modification of the Fei 
scheme 

Using the procedure in [9] one can derive the following symmetric 2-step scheme, 
i Yt + ( 2 U ™ 1 + (1 ~ 9)U " X + 2 U ^ 

_ ^7 \ TT n\2 ( JT n+l , TT n-l\ , M 1 ~ l) 

2 



\u n \ 2 (u n+1 + c/"- 1 ) + W—R (f/") 2 (u n+1 + u n - 1 ) 



9 and 7 are real parameters. Note that = 7=1 yields the method of Fei et 
al. This method is by construction energy preserving, and its energy function 
is given as 



m 

^E^I^T 



(l-7)fe 
16 

(l-7)fe 
16 



E A ((^ +1 ) 2 + ^™ +1 ) 2 ) (oc) 2 + (C 
E A ((^ +1 ) 2 - (C +1 ) 2 ) (to 2 - (O 2 ) 



The last two lines comes from the relation |u| 4 = | (it 2 + u 2 ) — | (u 2 — u 2 ) . We 
are only aware of a conserved density function in the Fei case, we will however 
see that some of the parameters yield improved stability compared to Fei. 
The dispersion relation is 

sin cut — K 2 (9 cosujt + (1 — 0)) + q cos cot. 

Using the same procedure as in chapter |2.2| we get the stability polynomial 

p(z) = ( Q c_, - (1 - 7 ) V) z 4 

+ (7(7 - 1)6 2 - ((2 - 7)6 + d t ) c-e - ((2 - 7)6 + d- t ) c t ) z 3 

+ (6(2 - j)(d e + d-t) + dtd-i + cic-e + c e c_ e + 2(1 - j) 2 q(b - q) + 4(1 - 7)6 

+ (7(7 - 1)6 2 - ((2 - 7)6 + dt) c-i - ((2 - 7)6 + d- t ) ct) z 

+ (c t c-e - (1 - 7) V) , 

where 

ct = (i - jq - 9{K + L) 2 ) e~ iuT and d e = 2(1 - 6){K + L) 2 . 

The scheme is stable provided that all zeroes of the self-reciprocal polynomial 
p(z) is on the unit circle. In figure 3.1 we see that the stability region is almost 



the same size as for Besse, compare with figure |2.4| These parameters clearly 
yield an improvement over the Fei method. 



10 



61 = 0.5, 7 = 1, K = 1 = 0.5,7 = 1 




Figure 3.1: Stability for the modified scheme. 



4 Conclusion 

We have seen that two schemes which have essentially the same conserving 
properties and are both symmetric and linearly implicit show a rather different 
behaviour when applied to the cubic Schrodinger equation. This behaviour is 
analysed in terms of plane wave solutions. Both schemes can be interpreted 
as two-step schemes and thus have a spurious solution, in the Fei scheme this 
solution is unstable even for small time steps contrary to the scheme of Besse. 
We show however that the scheme by Fei et al. can be stabilised without losing 
the symmetry or the energy preservation property. 

It is interesting to observe that two different schemes, designed to be con- 
servative in a very similar way, turns out behave completely differently with 
respect to stability. It remains to be seen if this phenomenon appears also in 
other PDE models than the cubic Schrodinger equation. 
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